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Abstract 

The domain size dependence of piezoelectric properties of ferroelectrics is investigated using a continuum 
Ginzburg-Landau model that incorporates the long-range elastic and electrostatic interactions. Microstruc- 
tures with desired domain sizes are created by quenching from the paraelectric phase by biasing the initial 
conditions. Three different two-dimensional microstructures with different sizes of the 90° domains are 
simulated. An electric field is applied along the polar as well as non-polar directions and the piezoelectric 
response is simulated as a function of domain size for both cases. The simulations show that the piezoelec- 
tric coefficients are enhanced by reducing the domain size, consistent with recent experimental results of 
Wada and Tsurumi (Brit. Ceram. Trans. 103, 93, 2004) on domain engineered BaTiO^ single crystals. 



1 



I. INTRODUCTION 



Ferroelectric s are excellent piezoelectric materials that can convert electrical energy into me- 
chanical energy and vice versa 1 . This electromechanical property arises due to the coupling of 
spontaneous polarization with lattice strain. Many devices such as ultrasonic transducers and 
piezoelectric actuators make use of this property 2 . Recently, there has been considerable interest 
in this field due to the observation of a giant piezoelectric response if the applied field is along a 
non-polar directional It is believed that this "superpiezoelectric" response is due to the symmetry 
change caused by a rotation of the polarization towards the direction of the applied field 5 . Domain 
configurations produced by the field applied in the non-polar direction are termed engineered do- 
mains. There are also a large number of domain walls between the degenerate variants which 
affect the piezoelectric property. In a recent paper, Wada and Tsurumi 6 studied the dependence of 
the piezoelectric properties of domain engineered BaTi0 3 single crystals as a function of domain 
size. Engineered domain configurations with a range of domain sizes were synthesized. The study 
revealed that piezoelectricity is enhanced for domain engineered crystals with small domain sizes 
(or high domain wall density). Thus, domain walls influence the piezoelectric properties and it is 
important to compute the contribution of the domain walls to the piezoelectric response. 

Electromechanical properties of ferroelectric s have been studied theoretically using first- 
principle calculations 5 -^. A continuum Landau theory describing a single domain or homoge- 
neous state has been used to study the electromechanical properties of BaTiO^ as a function of 
temperature and electric field direction 9 . Although such calculations provide valuable insights 
into the physics of the polarization- strain coupling, they do not describe inhomogeneities due to 
domains and domain walls. Recently, we studied the piezoelectric properties of domain engi- 
neered two-dimensional (2D) ferroelectric s using the time-dependent Ginzburg Landau (TDGL) 
theory 10 . The important conclusion from our simulations was the role played by the domain walls 
in nucleating an electric field induced structural transition if an electric field is applied along a 
non-polar direction. We showed that the field induced transition occurred at lower electric fields 
for a multi-domain state, compared to an analogous situation for a single domain state. To under- 
stand the recent experimental results of Wada and Tsurumi 6 that show piezoelectric enhancement 
at small domain sizes, we extend in this paper the TDGL model to investigate the dependence of 
piezoelectricity on the size of the 90° domains in the system and the domain wall density. Unlike 
Ref. [10] where the domain micro structure was obtained by quenching from the paraelectric phase 



2 



with random initial conditions, here we create domain structures with desired sizes by appropri- 
ately biasing the initial conditions. This procedure allows us to obtain domain microstructures 
with a range of domain sizes. The size dependence is studied for the case with the electric field 
along a polar axis as well as that with the field along a non-polar direction. 

The paper is organized as follows. In Sec. II, we describe the model in detail. Section III 
describes our simulations for the case of the electric field applied along a polar axis. In Sec. IV, 
we discuss the case in which the electric field is applied along a non-polar direction. We conclude 
in Sec. V with a summary and discussion. 

H. THE MODEL 

The calculations are based on a time-dependent Ginzburg-Landau mode l n i 12 i 13 with long-range 
elastic and electrostatic effects. We restrict ourselves to a 2D ferroelectric transition to illustrate 
the basic principles and use parameters from a model for BaTiO^ in our calculations 9 . The free- 
energy functional for a 2D ferroelectric system is written as F = Fi + F ern + F es . Here i 7 } is the 
local free energy 9 - that describes the ferroelectric transformation and is given by 



where P x and P y are the polarization components. The free energy coefficients a±, an, a U2 
determine the ferroelectric phase and the gradient coefficients gi, g 2 and g 3 are a measure of do- 
main wall energies. E x and E y are the components of an external electric field. Elastic properties 
are studied by using the strains r\\ = r/ xx + rj yy , r\ 2 = Vxx — Vyy an d Vs = Vxy, where 77^ is the 
linearized strain tensor defined as 77^ = (uij + Uj yi )/2 = x,y), Ui being the components 
of the displacement vector. The electromechanical coupling is described in terms of these strain 
variables with the free energy 



F em = X df [{77! - Q 1 (P X 2 + P y 2 )} 2 + {r/ 2 - Q 2 (P X 2 - P 2 )} 2 + {r/ 3 - Q 3 P X P V } 2 ] . (2) 



Here Q±, Q 2 and Q 3 are obtained from the electrostrictive constants of the material with Q\ = 
Qn + Q12, Q2 = Q11 — Q12 and Q 3 = Q44 (electrostrictive constants describe coupling be- 
tween strains and polarization, that is, rj xx = Q\\P 2 + Qi 2 P y 2 , r\ yy = Q\\P y + QnPx 2 and 
Vxy = QaPxPy)- Notice that the free energy F em vanishes for a homogeneous state since the 
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homogeneous strains in equilibrium are given by r]i e = Qi(P x 2 + P y 2 ), f] 2 e = Qi^Px' — Py 2 ), 
and r] 3 e = Q 3 P x P y . However, this free energy does not vanish for an inhomogeneous state. For an 
inhomogeneous state, the strains rji, r\ 2 and 773 are related to each other by the elastic compatibility 
constraint 

2 d 2 d 2 d 2 
^ ^dx 2 dy 2 ^ 2 dxdy 1 ^ 3 ^ ^ 
Using this relation, the strain 771 can be eliminated from F em resulting in a nonlocal interaction 

between the strains involving rj 2 and 773. Using the equilibrium strains defined by 7] 2 e and 7] 3 e , the 

electromechanical free energy can be written as 



F em = A / dk 



C 2 (k)T 2 (k) + C 3 (k)T 3 (k) -T^k) 
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(4) 



where the homogeneous state corresponding to the k = mode has been excluded from the above 
integral. The constant A is the strength of this nonlocal interaction and hence it influences the 
underlying microstructure. The quantities Ti(k), T 2 (k) and T 3 (k) are respectively the Fourier 
transforms of Q^P 2 + P 2 ) , Q 2 {P 2 - P 2 ) and Q 3 P x P y ; C 2 = {k 2 - k y 2 )/{k x 2 + k 2 ) and 
C 3 = k x kyl [k 2 + k y 2 ) axe the orientation dependent kernels. The electrostatic contribution to the 
free energy is calculated by considering the depolarization energy 15 

Fes = -nj df [E d -P + e (E d ■ E d /2)} , (5) 

where E d is the internal depolarization field due to the dipoles and \i is the strength of this interac- 
tion. The field E d can be calculated from an underlying potential using E d = — V0. If we assume 
that there is no free charge in the system, then V • D — 0, where D is the electric displacement 
vector defined by D = eoE d + P. This equation gives rise to the constraint — e o V 2 + V ■ P = 0. 
The potential <p is eliminated from the free energy F es using the above constraint to express F es in 
Fourier space as 



7 es = I 
1€q 



k x P x (k) + kyPy(Ji) 
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(6) 



The above integral excludes the homogeneous k = mode which means that the homogeneous 
depolarization field due to surface charges has been neglected. The total energy is defined by 
F = F t + F em + F es with two additional constants, i.e. A and /i are essential for the description of 
multi-domain states. 

The dynamics of the polarization fields is given by the relaxational time-dependent Ginzburg- 
Landau equations 



dP SF 



where 7 is a dissipation coefficient and i = x,y represents the polarization components. We first 
introduce rescaled variables defined with u = P x /Pq, v = P y /Po, C — r/5 and t* = 7|ai(T )|t, 
where T is a fixed temperature. In this work, we use the parameters^ for BaTi0 3 for the local 
part of the free energy F t . The parameters which can be dependent on the temperature T are: 
«i = 3.34 x 10 5 (T - 381) VmC" 1 , a n = 4.69 x 10 6 (T - 393) - 2.02 x 10 8 Vm 5 CT 3 , a in = 
-5.52 x 10 7 (T - 393) + 2.76 x 10 9 Vm 9 C~ 5 , a 12 = 3.23 x 10 8 Vm 5 C~ 3 and a U2 = 4.47 x 10 9 
Vm 9 C~ 5 . The electrostrictive constants are given as Qn = 0.11 m 4 C~ 2 , Q12 = —0.045 m 4 C~ 2 
and Q44 = 0.029 m 4 C~ 2 . We assume that the coefficients g 1 = g 2 = g% = g and use the value 
g = 0.025 x 10~ 7 Vm 3 /C quoted in the literature 16 . To calculate the rescaled quantities, we use 
T = 298K, Pq = 0.26 Cm~ 2 and S ~ 6.7 nm. The values chosen for the long-range parameters 
are A = 0.25|ai(T )|/P 2 and /1 = 20e |ai(T )|. All results presented below are expressed in 
terms of the rescaled time t*. 

HI. SIMULATIONS 

The time-dependent Ginzburg-Landau model with the above rescaled parameters is used to 
simulate the domain patterns and electromechanical properties. The equations are discretized on a 
128 x 128 grid with the Euler scheme using periodic boundary conditions. For the length rescaling 
factor 5 ~ 6.7 nm, this discretization corresponds to a system of size ~ 0.85 fim x0.85 /im. We 
simulate the properties of this 2D model at T = 298K. At this temperature, the minima of the 
free energy Fi define a rectangular ferroelectric phase with the four degenerate states (±0.26, 0) 
Cm -2 and (0, ±0.26) Cm -2 . Since we want to study the domain size dependence of properties, 
we create domain structures with required domain size instead of letting the domain structure form 
after a quench from the paraelectric phase. This is achieved by choosing initial conditions based 
on the following procedure. We consider a function 

R(x,y)= cos ms y (8) 

The initial conditions are set up by 

P x (x, y)=P , P y (x, y) = 0, R(x, y) > 

P x (x,y) = 0, P y (x,y) = P . R(x,y)<0 (9) 

These initial conditions ensure that multi-domain states with head to tail domain walls oriented 
along [11] are formed. The above initial conditions also ensure that only two of the four variants 
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with head to tail domain walls are formed in the multi-domain. The quantity N controls the 
number of domain walls and hence the domain size of the resulting microstructure. We consider 
the cases N = 2, 4, 10 corressponding respectively to 90° domain patterns with mean domain sizes 
L ~ 0.3/im, 0.15/xm, 0.06yum. The top left snapshots in Figs. 1, 2 and 3 represent the prepared 
zero field multi-domain states for L ~ 0.3fxm, 0.15/im and 0.06/im, respectively. These were 
obtained by solving Eqs. (7) for a time interval t* = 100 using the initial conditions given by Eqs. 
(8) and (9). A close look at the local dipoles within the domains shows that the polarization vectors 
are (P 0) A) and (A, P ), unlike the ideal single crystals which are described by (P , 0) or (0, P ). 
This means that within the domains, the polarization vectors are slightly rotated compared to the 
single crystals. As the domain size becomes smaller, the quantity A increases and the polarization 
vectors within the domains get increasingly rotated from the ideal [10] and [01] directions. This 
rotation is very strongly observed for the smallest domain size L ~ 0.06/im (Fig. 3(a)). This is 
due to the fact that the domain walls are closely spaced and the length of the diffuse interfaces is 
comparable to the domain width. 

To simulate the effect of an external electric field, the evolution equations are solved with a 
varying E. We consider two cases: (A) Field applied along the [01] direction, E = (0,E ), 
corresponding to the polar direction. (B) Field applied along the [11] direction, E = (^|, ;%), 
corresponding to a non-polar direction. 

A. FIELD APPLIED ALONG A POLAR DIRECTION 

We first study the traditional scenario when the electric field is applied along one of the polar 
directions. In the present simulations, we apply the field along the [01] direction which is a polar 
direction. The field is applied quasi-statically, i.e. in fixed increments of A.Et 01 ] = 0.92kV/ cm and 
we let the system relax for t* = 100 time steps after each change. Since [01] is a polar direction 
for the parameters used in the present simulations, the state (0, 0.26)Cm~ 2 is favored. Figures 1, 
2 and 3 show the electric field induced domain evolution for domain patterns with mean domain 
sizes L ~ 0.3/im, 0.15/xm and 0.06/im, respectively. It can be seen that domains aligned along 
the [10] direction switch towards the [01] direction thereby forming a single domain state for all 
the three cases. A comparison of the evolution in Figs. 1, 2 and 3 shows that a single domain 
state is established at smaller electric fields for domain patterns with a larger number of domain 
walls (or smaller domain sizes). For example, for the smallest domain size L ~ 0.06/im (Fig. 
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3), the single domain is established at an electric field E\ i] ~ 2.5 kV/cm (this electric field is 
much smaller than the electric field required to create single domain states for the domain patterns 
with L ~ 0.3/im, 0.15/im). Figures 4(a) and 4(b) show the variation of the average polarization 
with the applied field for the evolution depicted in Figs. 1, 2 and 3. For comparison, we also 
show the polarization vs. electric field response of a single domain state polarized along [01], 
i.e. P = (0,0.26) Cm" 2 , when an electric field is applied along the [01] direction. Figure 4(a) 
shows the evolution of the [10] component of the average polarization Pioi- At zero field, for the 
three cases shown in Figs. 1, 2 and 3, P[ 10 ] has a non-zero value due to coexisting domains and 
domain walls. Interestingly, this average value increases with decreasing the mean domain size. 
This increase can be attributed to the rotation of the dipoles within the domains. As discussed 
earlier, the polarization vectors within the domains are given by P = (P , A) or P = (A, P ). 
Since A increases with decreasing domain size, the average values over the multi-domain states 
also increase as the domain size becomes smaller. 

As the electric field is applied along [01], P[ W ] decreases to zero for all three cases due to 
the switching of domains polarized along [10] towards the [01] direction. As discussed earlier, 
Prio] reaches zero fastest for the smallest domain size, i.e. L ~ 0.06/im. The single domain 
P[io] remains zero, as expected. In Fig. 4(b), we plot the average polarization along the [01] 
direction, P[oi], as a function of the applied field E\n\\. Since [01] is a polar direction, P[ i] grows 
for all the cases. Here also, the case with the smallest domain size reaches the saturation value the 
fastest. P[ i] for the single domain varies only slightly with the applied field as there is no domain 
switching for that case. 

To study the electromechanical behavior, we have also computed the variation of the strains 
with the applied electric field. To evaluate the contribution of the applied electric field to the 
strain, we subtract off the zero field strain. Figures 5(a), 5(b) and 5(c) show the behavior of 
average strain components (rj xx (E [01] )} - (rj xx (E [01] = 0)), (r] yy (E [01] )} - (rj yy (E [01] = 0)) and 
(rj xy (Eioi]))-(r} xy (E[oi] = 0)) respectively. Here 77^. = Q n P x 2 +Q 12 P y 2 ,Vyy = QiiPy 2 +Qi2P x 2 
and 7] xy = Q^P x P y . Since the multi-domain states switch to a single domain state polarized along 
[01], shrinkage along the transverse [10] direction is observed, as can be seen in Fig. 5(a). The 
magnitude of this transverse strain is almost the same for all the multi-domain states, although 
the field required to establish the single domain state increases as the domain size is increased. 
The corresponding single domain undergoes very small shrinkage along the transverse direction 
as there is no domain switching involved. Figure 5(b) shows the behavior of the longitudinal 
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strain along the direction of the applied field for the three multi-domain states as well as the 
corresponding single domain states. An extension along the [01] direction is observed for all the 
cases. However, the multi-domain states generate much larger strains in comparison to the single 
domain state. This is due to the 90° domain switching in the multi-domain states that results in 
the extra strain. Figure 5(c) shows the average shear the crystal undergoes during the evolution 
depicted in Figs. 1, 2 and 3. It is clear that the magnitude of the shear depends on the number of 
domain walls in the system. This is due to the fact that the domain walls in the unpoled multi- 
domain states are sheared relative to the bulk. Upon applying the electric field, these domain walls 
disappear resulting in a net shear strain. Thus the average shear depends on the number of domain 
walls. The positive shear strain for low fields ("overshoot") is due to the domain switching process. 
Since there are no domain walls in the single domain state, shear strain is zero for all values of the 
electric field, as can be observed in Fig. 5(c). 

We have also studied the domain size dependence of the longitudinal piezoelectric coefficient 
c4g • The piezoelectric coefficients are calculated from the slope of (%/(-E|oi])) — (Vyy (E[oi] = 0)) 
vs. E[qi] curve in Fig. 5(b). Figure 6 shows the behavior of df^ vs. E[oi] for the three multi- 
domain cases along with the analogous single domain case. The high values observed in the 
electric field range — 10 kV/cm are due to the switching of domains. To clearly show the 
behavior of piezoelectric constants in the low-field regime, we replot the data of this figure in the 
inset for df^ < 1200. The data in the inset shows that the low-field piezoelectric coefficients 
are enhanced as the domain size is decreased. For example, for the smallest domain size Lq ~ 
0.06/xm, 4° 3 1] ~ 1100 pC/N compared to <ff ~ 210 for L ~ 0.3/zm. In the large field regime 
(E[ i] > 10 kV/cm), df^ is nearly equal for all the cases as they all correspond to a poled single 
domain state. 

B. FIELD APPLIED ALONG A NON-POLAR DIRECTION 

In this section, we study the case when the configurations depicted in Figs. 1(a), 2(a) and 3(a) 
are subjected to an electric field along the [11] direction. This situation is a 2D analog of the 
experiments by Wada and Tsurumi 6 where the electric field was applied along the [1 1 1] direction 
to tetragonal multi-domain single crystals of BaTiO^. 

Here, we apply a quasi-static electric field along the [11] direction. The field is applied in 
increments of AEVqi = 0.92 kV/cm and the configurations are allowed to relax for t* = 100 time 
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steps after each change. For this case, application of the field does not immediately result in the 
creation of a single domain state along [11]. Instead, the multi-domain structure remains stable 
and the polarization vectors rotate until an electric field induced transition to a [11] polarized 
rhombic state takes place. This situation is depicted in Figs. 7, 8 and 9 corresponding to the multi- 
domain states with domain sizes L ~ 0.3/im, 0.15/im and 0.06/xm, respectively. It is observed 
that the field induced transition occurs at a lower electric field as the domain size decreases. This 
result corroborates our earlier conclusion that the domain walls help nucleate the field induced 
transition 10 . Thus, the larger the number of domain walls, the smaller the field required to induce 
the transition. For example, for the smallest domain size L ~ 0.06/xm, the transition occurs 
at E[u] ~ 4.6kV/cm whereas for the largest domain size L ~ 0.3/im, the transition occurs at 
E[u] ~ 24 kV/cm, a more than 50% change. 

The evolution of the components of the average polarization for the situations depicted in Figs. 
7, 8 and 9 is plotted in Figs. 10(a) and 10(b). The response of a single domain state with initial 
polarization (0, 0.26)Cm~ 2 is also shown. The zero field components P[ 10 ] and Ppi] of the average 
polarizations for the multi-domain states in Figs. 7, 8 and 9 are non-zero due to the coexisting 
domains of (0.26, A)Cm~ 2 and (A, 0.26)Cm~ 2 . Since the polarization vectors rotate towards the 
[11] direction, the evolution of Pri ] and P i] is almost identical as both the polarization variants 
exist in nearly equal proportion. The single domain on the other hand starts from (0, 0.26)Cm~ 2 
till it transforms to a rhombic state (0.21, 0.21)Cm~ 2 . Figure 10 also shows that the field required 
to transform the multi-domain state to a rhombic phase depends on the number of domain walls in 
the system. However, the polarization components after the transition are same for all the cases as 
eventually a single domain rhombic state is established. 

Figures 11(a), 11(b) and 11(c) show the evolution of (^(-Ejii])) — {Vxx{E\\x\ = 0)), 
(r)yy(E[ n] )) - (r] yy (E [u] = 0)) and (77^(^11])} - (^(^[ii] = 0)), respectively. The results 
of this figure can be understood in terms of the electric field induced symmetry changes. Let us 
first examine the results for the single domain state. The zero field initial state (0, 0.26) Cm~ 2 
corresponds to a rectangular symmetry whereas the final state (0.21, 0.21)Cm -2 corresponds to a 
rhombic symmetry. This symmetry change is achieved by a uniaxial shrinkage along [01] and a 
uniaxial extension along [10], as can be inferred from Figs. 1 1(a) and 1 1(b) (notice the sharp jump 
near the field induced transition). The behavior of shear strain [shown in Fig. 1 1(c)] is governed 
by the rotation undergone by the polarization vector. In contrast, the zero field multi-domain states 
correspond to a nearly square macroscopic symmetry due to the coexistence of two polarization 
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(rectangular) variants. Hence, the multi-domain evolutions of Figs. 7, 8 and 9 effectively corre- 
spond to electric field induced square to rhombic transitions. The jump due to the field induced 
transition occurs at smaller electric field values as the domain size is decreased. The uniaxial strain 
the crystal undergoes after the transition is almost the same along the [10] and [01] directions. In- 
terestingly, the saturation value of the strains is essentially the same for all the three multi-domain 
evolutions. The magnitude of shear strains after the transition, on the other hand, depends on the 
domain size (or the number of domain walls) in the initial state. As seen in Fig. 1 1(c), the amount 
of shear experienced by the crystal is the largest for the case with L ~ 0.3/xm and the smallest 
for the Lq ~ 0.06/im case. Pre-existing shear strains at the domain walls limit the total shear 
experienced by the multi-domain crystals and thus the larger the number of domain walls in the 
initial state, the smaller the shear strains produced. 

Figure 12 depicts the behavior of the longitudinal piezoelectric coefficients for the three 
multi-domain states as well as the analogous single domain situation. The quantity is calcu- 
lated from the slope of the longitudinal strain resolved along the [11] direction vs. Eny curve. 
The resolved strain is calculated as (r)[n](E[n])) -{r)[n] (E[n] = 0)), where r)[n] is given by 

run] = ^(Vxx + Vyy + Vxy)- (10) 

It is clear that the low-field piezoelectric coefficients for the smallest domain size L ~ 0.06/xm 
are more enhanced compared to the domain patterns with L ~ 0.15/xm and 0.3/im. The low 
field piezoelectric coefficients for the coarser domain patterns are not much higher than the single 
domain coefficients, consistent with recent experiments 17 . We believe that the enhancement is 
related to the response of unit cells in the domain wall regions as such regions become bigger as 
the domain size becomes smaller. 



IV. SUMMARY AND DISCUSSIONS 



We have used a Ginzburg-Landau formalism to study the domain size dependence of the piezo- 
electric properties. The present work is inspired by the recent experiments of Wada and Tsurumi 6 
on domain engineered BaTiO^ single crystals where the effect of the size of non-180 domains 
on the piezoelectric constants was studied. In our model calculation, we solved the 2D time- 
dependent-Ginzburg-Landau equations^ with biased initial conditions (the free energy parameters 
for BaTiOz were chosen from Ref. 0) to create three different multi-domain states with different 
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domain widths. 

Two different directions of the applied field were considered. In the first case, the multi-domain 
states were subjected to an electric field along the [01] direction, which is one of the four polar 
directions. The multi-domain states switched to single domain states polarized along the [01] 
direction, with the state having the largest number of domain walls switching at the lowest electric 
field. The multi-domain state with the smallest domain size also exhibited the largest value of the 
longitudinal piezoelectric constant df^ . This enhancement of the piezoelectric coefficient as the 
domain size is decreased reflects the metastability of the multi-domain states which can become 
easily switchable as the number of domain walls is increased. However, this enhancement of the 
piezoelectric constant may not be very useful in practical applications as the multi-domain states 
are not stable over a large range of electric fields. 

We also considered the case where an electric field along the [11] direction is applied to the 
same multi-domain states. This situation is analogous to the experiments of Wada and Tsurumi 6 
who studied the piezoelectric properties of tetragonal domain engineered BaTi0 3 single crystals 
under an electric field along [111]. For this case, we found that the multi-domain states remain 
stable until a field induced rectangular to rhombic transition takes place. Interestingly, we found 
that the transition occurs at smaller fields as the domain size is decreased. Since proximity to the 
field induced transition enhances the piezoelectric constants, the low-field piezoelectric constant 
for the smallest size simulated by us is found to be significantly higher than that for the single 
crystal and multi-domains with bigger domain sizes. Thus, the role of the domain walls in nucle- 
ating a field induced transition may be the cause of the enhanced piezoelectricity in small sized 
engineered domains observed by Wada and Tsurumi 6 . This enhancement may be used in practical 
applications, provided the field is not too close to the field induced transition. 
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Figure Captions: 

Figure 1 : Evolution of domains for an electric field applied along the [01] direction to the multi- 
domain state with domain size L ~ 0.3[im. The corresponding electric field levels are indicated 
at the top of each snapshot. 

Figure 2: Evolution of domains for an electric field applied along the [01] direction to the 
multi-domain state with domain size L ~ 0.15/zm. The corresponding electric field levels are 
indicated at the top of each snapshot. 

Figure 3: Evolution of domains for an electric field applied along the [01] direction to the 
multi-domain state with domain size L ~ 0.06/xm. The corresponding electric field levels are 
indicated at the top of each snapshot. 

Figure 4: Evolution of average polarizations P[i ] (Fig. 4a) and P[ i] (Fig- 4b) with the applied 
field E[oi]. The lines with circles correspond to the multi-domain state of Fig. 1, lines with crosses 
correspond to the multi-domain state of Fig. 2 and lines with squares correspond to multi-domain 
state of Fig. 3. Solid lines correspond to the single domain state. 

Figure 5: Evolution of average strains (rj xx (E [01] )} - (r]xx(E [01] = 0)) (Fig. 5a), (r] yy (E [01] )) - 
(rjyy(E [01] = 0)) (Fig. 5b) and (r) xy (E [01] )) - (rj xy (E [01] = 0)) (Fig. 5c) with the applied field 
E\ i]. The lines with circles correspond to the multi-domain state of Fig. 1, lines with crosses 
correspond to the multi-domain state of Fig. 2 and lines with squares correspond to multi-domain 
state of Fig. 3. Solid lines correspond to the single domain state. 

Figure 6: Variation of df^ (the longitudinal piezoelectric constant along [01]) with E\ i]. The 
lines with circles correspond to the multi-domain state of Fig. 1, lines with crosses correspond to 
the multi-domain state of Fig. 2 and lines with squares correspond to multi-domain state of Fig. 3. 
Solid lines correspond to the single domain state. The inset plots the data in the range df^ < 1200 
pC/N to show the low-field behavior. 

Figure 7: Evolution of domains for an electric field applied along the [11] direction to the multi- 
domain state with domain size L ~ 0.3/im. The corresponding electric field levels are indicated 
at the top of each snapshot. 

Figure 8: Evolution of domains for an electric field applied along the [11] direction to the 
multi-domain state with domain size L ~ 0.15/zm. The corresponding electric field levels are 
indicated at the top of each snapshot. 
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Figure 9: Evolution of domains for an electric field applied along the [11] direction to the 
multi-domain state with domain size L ~ 0.06/um. The corresponding electric field levels are 
indicated at the top of each snapshot. 

Figure 10: Evolution of average polarizations Pr 10 i (Fig. 10a) and P^i] (Fig- 10b) with the 
applied field Pmi . The lines with circles correspond to the multi-domain state of Fig. 7, lines with 
crosses correspond to the multi-domain state of Fig. 8 and lines with squares correspond to the 
multi-domain state of Fig. 9. Solid lines correspond to the single domain state. 

Figure 11: Evolution of average strains (r]xx(E [u] ))-(rj xx (E [n] = 0)) (Fig. 11a), (^(^u]))- 
(f]yy{E [u] = 0)) (Fig. lib) and (^(^[li])) - {r) xy {E [n ] = 0)) (Fig. 11c) with the applied field 
Prn]. The lines with circles correspond to the multi-domain state of Fig. 7, lines with crosses 
correspond to the multi-domain state of Fig. 8 and lines with squares correspond to the multi- 
domain state of Fig. 9. Solid lines correspond to the single domain state. 

Figure 12: Variation of df^ (the longitudinal piezoelectric constant along [11]) with Smi. The 
lines with circles correspond to the multi-domain state of Fig. 7, lines with crosses correspond to 
the multi-domain state of Fig. 8 and lines with squares correspond to the multi-domain state of 
Fig. 9. Solid lines correspond to the single domain state. 
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